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SECTION  I 


INTRODUCTION 


In  this  report  we  consider  the  problem  of  computing  sound  propagation 
in  ducts.  These  ducts  are  models  of  the  interior  of  jet  engines.  Since 
the  largest  source  of  jet  engine  noise  is  fan  noise,  there  has  been  an 
increase  In  interest  in  duct  acoustics  in  recent  years.  We  will  consider 
ducts  of  uniform  cross  section  with  or  without  mean  flow  and  ducts  with 
variable  cross  section  not  containing  flow.  The  ducts  are  either  two 
dimensional  or  axially  symmetric  three  dimensional.  To  compute  the 
sound  field  within  one  of  these  ducts,  the  momentum  and  continuity  equa- 
tions are  linearized  and  combined  to  yield  a linear,  second  order  partial 
differential  equation  in  the  acoustic  pressure.  The  boundary  conditions 
consist  of  a specified  pressure  distribution  at  the  duct  entrance,  a 
vanishing  normal  derivative  at  the  centerline,  and  a specified  acoustic 
impedance  at  the  exit  and  at  the  outer  wall  of  the  duct  which  corresponds 
to  an  acoustic  lining.  These  equations  are  then  solved  using  the  finite 
difference  method.  The  large  system  of  algebraic  equations  which  result 
from  using  finite  differences  is  solved  using  an  out-of-core  block  tri- 
diagonal system  solver  on  the  CDC  Cyber  74  computer.  The  acoustic 
velocities  are  then  determined  by  Integrating  the  linearized  momentum 
equations  and  then  used  to  compute  the  acoustic  power  at  the  entrance  and 
exit  of  the  duct.  These  are  then  used  to  determine  the  attenuation 
resulting  from  a particular  wall  lining. 

For  a nonuniform  duct  the  additional  step  of  mapping  the  nonuniform 
geometry  to  a uniform  duct  is  required  before  the  finite  difference 
method  is  applied  to  obtain  the  solution.  We  use  a conformal  map  so 
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that  right  angles  in  the  original  coordinate  system  become  right 
angles  in  the  mapped  coordinate  system.  The  mapped  equations  are  then 
solved  using  finite  differences  as  indicated. 

In  Section  II,  we  will  describe  in  more  detail  the  equations 
and  boundary  conditions  which  need  to  be  solved  and  give  the 
expressions  used  to  compute  the  acoustic  power.  In  Section  III, 
we  discuss  the  conformal  mapping  procedure  and  derive  the  mapped 
equations  and  boundary  conditions.  In  Section  IV,  we  derive  the 
finite  difference  equations  for  both  the  original  and  mapped  equations, 
and  indicate  the  method  used  to  solve  the  resultant  system  of  algebraic 
equations.  In  Section  V,  we  describe  the  computer  program  DUCT 
including  minor  program  modifications  for  extra  capability,  and  in 
Section  VI  we  give  samples  with  input  and  output. 

SECTION  II 
ANALYTIC  MODEL 


Differential  Equations 

If  viscous  and  heat  transfer  terms  are  neglected,  the  equations 
of  motion  for  an  ideal  gas  are 


Momentum 

p'  D(V')/Dt  = -Vp' 

(la) 

Continuity 

D(p’)/Dt  + p’(V*V')  = 0 

(lb) 

Energy  and  state 

p'  = const  p'*^ 

(Ic) 

where  p'  is  the  density,  V is  the  velocity  vector,  and  p'  is  the 
pressure.  If  it  is  assumed  that  the  pressure  perturbations  are 
small  compared  to  the  average  pressure  and  that  the  steady  velocities 
are  negligible,  Eq.  (1)  becomes 
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(2a) 


ap/3t  = -p^c^V‘V 

3V/8t  = (2b) 

When  combined,  Eq$.  (2)  reduce  to  the  wave  equation 

= (l/C^)  92-/3t2  (3j 

and,  with  time  dependence  of  the  form  [that  is,  p(X,Y,Z,t) 

= e^^^p(X,Y,Z)],  Eq.  (3)  becomes  the  Helmholtz  equation  for  the 

pressure  p 

+ (a)^/C^)p  = 0 (4) 

Analyzing  only  one  frequency  for  fan  noise  is  reasonable  since 

noise  at  different  frequencies  can  be  superimposed.  In  nondimensional 

12  3 4 

coordinates,  » » » Eq.  (4)  becomes 


7%  + (2Tm)^P  =0  (5) 

with  n = Ru)/(27rc)  for  a cylindrical  duct  (R  being  the  duct  radius), 
f)  = Ha)/{2TTc)  for  a rectangular  duct  (H  being  the  duct  height).  In 
two  dimensions,  for  rectangular  ducts,  Eq.  (5)  is 

3^p/3x^  + 3^p/3y^  + ( 27171  )^p  = 0 (6) 

whereas,  in  three  dimensions  with  axial  symmetry,  Eq.  (5)  becomes 

3S/3x^  + 3^p/3r^  + (l/r)  3p/3r  + (27Tn)S  = 0 (7) 

where  x is  the  axial  direction  and  r is  the  radial  direction.  In 
the  mean  flow  case,  the  steady  velocity  in  the  axial  direction  is 

assumed  to  be  a constant  U (and  no  longer  negligible),  yielding  the 

3 

dimensionless  axially  symmetric  reduced  wave  equation 
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(8) 


3^p/3x^  + 3^p/3r^  + (l/r)  3p/9r 

- 3p/3x  + (21^1  )^p  - 0 

where  M = U/c  is  the  Mach  number  and  n * Rw/(2ttc)  is  the  dimensionless 
frequency,  R being  the  duct  radius. 

Boundary  Conditions 

To  absorb  acoustic  energy  within  the  duct,  the  outer  walls  are 
lined  with  Helmholtz  resonators.  The  Helmholtz  resonator  can  consist 
of  a perforated  sheet  placed  at  distance  d from  a solid  backing  sheet. 
The  acoustic  impedance  of  this  wall  lining  is  determined  by  the  hole 
diameters  and  backing  depth. 

If  Z is  the  acoustic  impedance  at  the  outer  wall,  then 

Z = p/n.V 

where  n is  the  exterior  normal  to  the  outer  wall.  For  no  flow, 
iwp^V  = -Vp,  so 

iZ/top^  = p/n*Vp 

(9) 

or 

3p/3n  + iwpjj  p/Z  = 0 

at  the  wall.  [Note  that  3p/3n  = n-Vp  by  definition.]  For  the  case 
of  uniform  flow  the  continuity  of  particle  displacement  is  used  rather 
than  the  continuity  of  particle  velocity  to  derive  the  dimensionless 
boundary  condition 

3p/3n  -2Tmip/Z  - (2M/Z  ) 3p/3x  + iM^/(27rTi7  ) 3^p/3x^  {10} 
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at  the  outer  wall  (see  Ref.  5 for  the  derivation).  For  a variable  wall 
lining  {such  as  a three-sectional  lining)  is  a function  of  x instead 
of  being  constant,  as  is  the  case  for  a uniform  lining.  The  remaining 
boundary  conditions  for  a cylindrical  duct  are 


p(0,r)  = f(r) 

at  the  entrance 

(lla) 

3p/3n  = 0 

at  r = 0 (the  centerline) 

(lib) 

3p/9n  + Zirnip/Z  = 0 

at  the  exit 

(11c) 

For  the  rectangular  duct  the  only  change  in  the  boundary  condition 
is  that  p{0,y)  = f(y)  is  specified  as  an  entrance  condition. 

Sound- Power  Attenuation 

For  both  the  no-flow  and  uniform  flow  cases,  the  expression  used 

3 4 

to  compute  the  dimensionless  axial  intensity  ’ is 

l(x,r)  = real  {p*u)/(27rTi)  + M[p*p  + (u*u  + v*T)/(2Trn)  ]/2  p2) 

where  u is  the  velocity  in  the  axial  direction,  v is  the  velocity 
in  the  radial  direction,  and  the  asterisk  denotes  complex  conjugation. 
For  the  no  flow  case,  of  course,  M = 0 and 


I = real  (p*u)/(2Trn) 


(13) 


To  compute  u and  v,  once  p is  known,  the  dimensionless  x-momentum  and 
r-momentum  equations  are  used: 


u = 1 3p/3x  + i M/(27rri)  3u/3x  (14a) 

V = 1 ap/3y  +■  i M/(27m)  3v/3x  (14b) 
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The  initial  condition  used  to  solve  for  u is  that  u = Zttop  at  the 
exit.  Then  the  irrotational ity  condition  3u/3y  = 3v/9x  and 
Eqs.  (14)  are  used  to  solve  for  v.  Once  the  intensity  is  known 
the  dimensionless  acoustic  power  is  computed  as 


E 

X 


l(x,r )rdr 


for  a cylindrical  duct,  or 


E = f I(x,y)dy 

Jo 


(15) 


(16) 


for  a rectangular  duct, 
as 


The  sound  attenuation  is  then  determined 


as  = 10  log^j,  (E^/EJ 


(17) 


where  E is  the  sound  power  at  the  entrance  and  E is  the  sound 

0 X 

power  at  the  axial  point  x. 


SECTION  III 
CONFORMAL  MAPPING 


The  well-known  Riemann  Mapping  Theorem  establishes  that  an 
arbitrary  simply  connected  domain  in  the  plane  can  be  mapped 
conformally  onto  an  open  rectangle  with  vertices  (0,1),  (0,-1), 
(X,-l),  and  (X,l).  (For  a diagram  of  such  a map,  see  Fig.  1.  Other 
examples  of  variables  of  ducts  are  given  in  Figure  2.)  If  the  map 
is  of  the  form 

z = z(x,y);  r = r(x,y) 

then,  for  it  to  be  conformal,  the  Cauchy-Riemann  equations 
3z/3x  = 3r/3y;  3z/3y  = - 3r/3x 


must  be  satisfied. 
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Assuming  axial  symmetry,  Eq.  (4)  may  be  rewritten  as 
3^p/3z^  + 3^p/3r^  + (l/r)  3p/3r  + (w/c)^p  = 0 


where  z is  the  axial  coordinate  and  r is  the  radial  coordinate. 
In  transformed  coordinates  this  equation  becomes 


0 = 


1 32Cx,y3  3g 

rCx,yJ  iy  3x 


+ 


1 

rCx.y) 


3zCx,y)  ^ 
3x  3y 


+ C2Trn3^ 


'3zCx,y)^ 

9x 


2 


+ 


(18) 


where  p(x,y)  = p(z,r)  and  n = Ra)/(27rc)  is  the  dimensionless  frequency. 
The  boundary  conditions  (9)  and  (11)  become 


]^C0,y)  = £(2(0,/)  ,rC0,y))  = gCy)  at  the  entrance  (19a) 


(19b) 

(19c) 

(19d) 


Finite  difference  solutions  will  be  obtained  in  Sec.  IV.  (See  Ref.  6 
for  a more  complete  description.) 


Exampl e 

For  positive  real  a < j,  the  cone  enclosed  in  the  (z,r)  plane  by 
the  lines  r = (sin  a)z/cos  a and  r = -(sin  a)z/cos  a and  by  arcs  of 
the  circles  z^  + r^  = and  z^  + r^  = R^e^"'’'  is  the  region  considered. 
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The  transformation  defined  by 


, o ax f,  ax  . 

z = Re  cos  ay;  r = Re  sin  ay 


(20) 


maps  the  open  rectangle  onto  the  prescribed  cone.  With  z(x,y) 
and  r(x,y)  defined  by  Eq.  (20),  the  finite  difference  equations 
corresponding  to  Eqs.  (18)  and  (19)  were  solved  on  a rectangle  using 
the  method  described  above. 

SECTION  IV 

FINITE  DIFFERENCE  SOLUTIONS  OF  ANALYTIC  MODEL 

To  obtain  solutions  of  Eqs.  (6-8)  in  the  case  of  a rectangular 
duct,  a no-flow  cylindrical  duct,  or  a cylindrical  duct  with  uniform 
flow,  respectively,  the  method  of  finite  differences  was  used.  The 
boundary  conditions  (10)  and  (11)  must  be  included  for  the  problems 
to  be  well  posed.  The  finite  difference  equations  may  be  summarized 
as  follows  (see  Refs.  2,  3,  6,  and  especially  7).  Set 


6x  = X/(n+l ) 
6r  = l/(m+l ) 


*0  ' 

r,,  = k6r 


i = 1 


n 


k = 1 , • • • , ra 


and 


where  m and  n are  positive  integers.  Then  for  a cylindrical  duct 
the  difference  equation  at  a point  not  adjacent  to  the  boundary 


* *®J,k+l  ■ ^■’j.k  * * *^J,k+l  " Pj,k-l'^*^*^’^'’'k' 

- - Pj.i,k>/<2''*’  " (2^>\.k  = 


Collecting  terms,  the  equation  becomes 

[2  + 2(l-M^)(6r)^/(6x)^  - (2TTn6r)^]  - Pj  l/(2k)] 

“ k-1^^  ■ “ 2TrnMi6r^/6x] 

■-  p [l-M^  + 2irnMi6r^/6x]  = 0 for  1 < J < n-1,  2 <_  k m-1. 


Similar  expressions  are  obtained  for  j=n,  k=l,  k=m,  which 
correspond  to  the  boundary  conditions  (10)  and  (11). 

To  obtain  a solution  of  Eqs.  (18)  and  (19)  on  a cylinder,  the 
following  finite  difference  scheme  is  employed. 

Set 

6x  = X/Cn  + 1), 

6y  = 1/  (m  + 1)  , 

Xj  = j 6x,  j = 1,  , n, 

= k 6y,  k = 1,  • • • , m, 


= If 

"j  ,k  " S ’>'P  ■ 

•■j  .k  ■ ’’'P  ’ 

and 


j = 1,  n, 
k = 1,  *•*,  m, 


Pj  ,k  ■ ’>'P  ’ 
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where  m and  n are  positive  integers.  Then  the  difference  equation 
is 


Collecting  terms,  the  equation  becomes 


’=  0 for  1 < j < n - 1,  2 < k < m - 1. 


Similar  expressions  are  obtained  for  j = n,  k = 1,  k = m,  which 
correspond  to  the  boundary  conditions  (19). 

For  a reasonable  number  of  mesh  points,  say  50  in  the  x 
direction  and  20  in  the  y direction,  the  linear  system  of  equations 
to  be  solved  consists  of  1000  unknown  p.  . *s  to  be  determined.  More- 

J , K 

over,  the  p.  . 's  are  complex  numbers  because  of  the  boundary 
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conditions,  so  there  are  2000  unknown  real  numbers  to  compute. 
Standard  iterative  methods  will  not  work  for  this  problem,  for  the 
system  is  neither  positive  definite  nor  symmetric.  Therefore,  a 
direct  method  was  used  to  solve  the  system,  taking  advantage  of 
its  sparseness. 

The  set  of  difference  equations  forms  a system  of  linear 
equations  which  can  be  written  in  matrix  form  as  MP  = F 

where  M is  the  finite  difference  matrix,  P is  the  unknown  pressure 
matrix  and  F is  the  matrix  containing  the  prescribed  pressures. 

M has  the  following  form: 

0...0-a0000  ...0-1  b-1  0000  ...  O-cOOOO  ...  0 
0 ...  0 0-a  0 0 0...  0 0 ^ bJOOO...O  0 -c  0 0 0 . . . 0 

0...0  00-a00,..0  00  -lb-1  00...  0 00-C00...0 

0...0  000-a0...0  00  0 -lb-1  0...0  000-C0...0 

which  can  be  written  as  a block- tridiagonal  matrix 
0 0 ...  0 

B2  0 ...  0 where  A,  B,  C are  square  matrices 

° ^3  ^3  ° which  fit  into  the  pattern. 

0 ... 

M can  now  be  factored  as 


M = L * U 


or 


I 0 ...  0 


0 0 ...  0 


0 ...  I 


0 U2  0 ...  0 


0 ...  0 


2 
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where 


U,  = B 


for  m = 2 , ...  n 


Therefore,  the  system  to  be  solved  may  be  written  as 


L * U * P = F 


Solving  L * Y = F for  Y and  then  solving  U * P = Y for  P gives 
the  desired  pressures. 


SECTION  V 


PROGRAM  DESCRIPTIONS 


Program  DUCT  computes  the  pressure  within  cylindrical,  rec- 
tangular, and  nonuniform  ducts  using  the  finite  difference  method 
described  in  Section  IV.  Using  these  values  for  the  pressure,  DUCT 
then  calculates  the  acoustic  power  and  attenuation. 

The  logical  flow  of  DUCT  is  represented  schematically  in 
Figure  3.  A summary  of  the  various  subroutines  called  by  DUCT  is 
now  given. 

DUCT  The  main  routine  which  controls  the  reading  of  input 


and  the  initialization  of  variables,  DUCT  calls 


FUNCT  directly  or,  if  the  function  is  to  be 


minimized,  calls  an  optimization  routine  ZXMIN. 


FUNCT 


A statement  function  which  does  the  calling  to  the 
subroutines  which  set  up  the  system  to  be  solved. 
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ZXMIN 


SETUP! 


SETUP 


solve  the  system,  compute  the  attenuation  of  sound, 
and  print  the  output.  By  having  a statement  function 
do  the  calling,  standard  optimization  routines  can 
be  used  which  minimize  a given  function. 

A library  subroutine  which  finds  the  minimum  of  a 
function . 


This  subroutine  uses  a finite  difference  method  to 
transform  the  wave  equation  for  a rectangular  duct 

P + P + * 0 

XX  yy 

into  centered  difference  equations  at  each  point. 
For  a point  not  on  the  boundary,  the  equation  has 
the  following  form: 


'■’’i+l.j 


Collecting  terms,  the  equation  becomes 
2 

^ 2 


P,.j- 


^ P.  , . - ^ P..,  . = 0 
j2  1+1  .J  0,2  1+1  .J 

X 


This  system  of  equations  is  then  written  in  matrix 
form. 


Using  the  methods  and  equations  presented  in  Section  IV, 
SETUP  produces  the  finite  difference  for  a nonuniform 
duct. 
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SETUP2 

This  subroutine  sets  up  the  matrix  for  a cylindrical 

duct  as  described  in  Section  IV. 

BLOCK 

To  increase  the  efficiency  of  the  program,  the  matrix 

produced  in  the  SETUP  routines  is  stored  in  block- 

tri diagonal  matrix  form.  BLOCK  calls  LU3  to  factor 

the  matrix  into  bi diagonal  form.  LEQS3  is  then 

called  to  solve  the  system  using  back  substitution. 

BLOCK  also  invokes  CLOSMS,  OPENMS,  WRITMS  and  READMS 

which  are  system  subroutines  available  on  Control  Data 

machines. 

OUTPUT 

Prints  the  real  and  imaginary  parts  of  the  pressure 

calculated  in  the  BLOCK  subroutine. 

SOUND 

Given  the  pressure,  SOUND  computes  and  prints  the 

sound  power  and  attenuation 

dB  = 10  log^Q  (E2/E^) 

where  E-j  is  the  acoustic  power  at  the  entrance  of 
the  duct  and  E2  is  the  acoustic  power  at  the  exit. 

Input 

The  information  needed  to  perform  the  calculations  can  be  entered 
through  cards  or  disks.  A description  of  required  input  is  presented 
here. 


Card  No. 

Columns  Variable  Name  Description 

1 

1-3  I DUCT  An  integer  representing  the 

duct  geometry.  Enter  -1 
for  a nonuniform  duct,  0 for 
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Card  No.  Columns  Variable  Name 


2 1-10  XM 


3 1-20  ZETAX 


4 1-20  ZETAX 


5 

1-10 

XMAX 

5 

11-20 

YMAX 

6 

1-10 

ETA 

7 

1-3 

I END 

Description 

a rectangular  duct  or  1 for 
a cylindrical  duct.  (13) 

Mach  Number.  For  the  acoustic 
linearization  to  be  valid 
-.5  £ XM  _<  .5.  Set  XM  equal  to 
zero  for  no  flow  case.  (FlO.O) 

Complex  number  representing  the 
exit  impedance  (FlO.O,  FlO.O). 

Wall  impedance  at  Y = YMAX, 

ZETAX  = 9999999999  + 19999999999 
corresponds  to  a hard  wall  (FlO.O, 
FlO.O). 

Length,  L,  of  duct  in  meters 
(FlO.O). 

Height,  H,  of  duct  in  meters.  If 
problem  has  been  made  non-dimensional, 
set  YMAX  = 1 and  XMAX  = (FlO.O) 

Frequency  parameter  (FlO.O). 

Number  of  pts  in  x-direction, 
excluding  pts  on  the  boundaries. 
Maximum  of  100  if  subroutines  READMS, 
WRITMS  are  used.  If  XM  0,  TEND 
should  be  odd.  (13) 
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Card  No. 

7 

8 

Set  KOUT  = 0 to  suppress  the 
printing  of  pressures  and 
KOUT  = 1 for  complete  output. 
(13) 


Col umns 

Variable  Name 

Description 

4-6 

I END 

Number  of  pts  in  y-direction. 

excluding  pts  on  boundaries. 

Maximum  of  20.  (13) 

1-3 

KOUT 

Controls  the  form  of  the  output 

If  ZXMIN  is  to  be  called  to  find  the  minimizing  values  of  ZETAY,  the 
following  information  is  required. 


Card  No.  Columns  Variable  Name 


Description 


4-6 


NUMMAX  Maximum  number  of  iterations  or 

calls  to  FUNCT  by  ZXMIN.  (13) 


7-16 


IH 


Number  of  digits  of  accuracy 
desired  in  the  estimates  of  ZETAY.  (13) 


If  KOUT  = 0,  the  following  card  can  be  left  blank. 


Card  No. 


Columns 


1-3 


Variable  Name 


lOUT 


Description 

Increment  value.  The  pressure 
will  be  printed  for  pts  with 
abscissas  corresponding  to  the 
values  of  I = 1 through  I END  in 
steps  of  lOUT.  (13) 
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Card  No. 


Columns  Variable  Name 


9 4-6  JOUT 


Description 

Increment  value.  The  program 
will  print  the  pressure  at  pts 
located  on  grid  lines  corresponding 
to  J = 1 through  JEND  in  steps 
of  JOUT.  (13) 


9 7-9  IT  The  pressure  will  be  printed  for 

pts  located  on  grid  lines  corre- 
sponding to  J = 1 through  IT. 
Maximum  value  is  IT  = JEND. 
(Headings  will  be  distorted  if 
JOUT  =1  and  IT  > 15). 


10  1-10  ALPHA  Generating  angle  in  radians  of  a 

nonuniform  duct.  Required  only 
if  IDUCT  = -1  . (FI 0.0) 


Multi  sectional  Linings 

If  the  user  desires  a multisectional  lining,  several  changes  must 
be  made  to  the  program  to  accommodate  the  additional  impedance  parameters 
(See  Section  IV).  The  following  changes  were  used  to  modify  the  program 
for  a three-sectional  lining  in  a nonuniform  duct.  Modifications  for 
rectangular  and  cylindrical  ducts  are  enclosed  in  brackets. 

Recall  that  the  impedance  is  a complex  number  so  that  an  n-sectional 
lining  means  ZET,  the  array  containing  these  values,  will  be  of  dimension 
2n.  To  adjust  the  dimension  of  ZET,  insert  the  appropriate  value  in 
line  7 of  DUCT.  In  DUCT,  arrays  WA,  HH,  and  G are  parameters  required 
by  ZXMIN  and  their  dimension  depends  on  the  dimension  of  ZET.  Set  the 


17 


dimension  of  WA  equal  to  N(N+4).  The  dimension  of  HH  equal  to  N{N+1)/2, 
and  the  dimension  of  G equal  to  N where  N is  the  dimension  of  ZET 
(N=2n).  For  example,  for  a 3 section  lining,  DUCT. 7 should  now  read: 

DIMENSION  ZET(6),  WA(18),  HH(21),  G(6) 

In  addition,  a dimension  statement  will  need  to  be  inserted  after 
line  8 in  SETUP  [SETUPl , SETUP2]  which  reads  DIMENSION  ZET  (6).  Also, 
in  line  55  of  DUCT,  the  number  following  FUNCT  in  the  parameter  list 
is  the  dimension  of  ZET,  thus  it  should  be  changed  to 

ZXMIN  (FUNCT,  6,  IH,  NUMMAX,  ZET,  HH,  G,  FMIN,  WA,  lER) 

To  initialize  the  values  of  ZET  to  those  of  a uniform  lining,  the 
following  cards  should  be  inserted  after  line  51  in  DUCT 

ZET(3)  = ZET{1) 

ZET(4}  = ZET{2) 

ZET(5)  = ZET{1) 

ZET(6)  = ZET{2) 

If  the  program  is  restarted,  the  expressions  on  the  right  side  of 
the  equations  can  be  replaced  with  the  minimizing  values  from  the 
previous  run. 

In  order  to  pass  the  additional  values  to  the  subroutines,  lines 
10  through  12  in  FUNCT  should  be  replaced  with 

IF(IDUCT.LT.O)  CALL  SETUP  (ZET) 

IF(IDUCT.EQ.O)  CALL  SETUPl (ZET) 

IF(IDUCT.GT.O)  CALL  SETUP2(ZET) 

and  the  first  line  of  SETUP  [SETUPl,  SETUP2]  should  be  replaced  with 
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SUBROUTINE  SETUP  (ZET) 

SUBROUTINE  SETUP! (ZET) 

_SUBROUTINE  SETUP2(ZET)_ 

To  print  the  values  of  ZET  before  each  pass,  insert  these  cards  after 
line  9 in  FUNCT 

PRINT  999,  (ZET{I),  I = 1 .6) 

999  FORMAT  (IX,  6F16.8) 

where  I ranges  from  1 to  the  dimension  of  ZET. 

Finally,  the  location  of  the  various  linings  are  inserted  after 
SETUP  .34  [SETUP1.22,  SETUP2.22].  For  example,  if  each  lining  occupies 
one  third  of  the  duct  the  following  lines  are  added; 

ZETAY  = ZET(1)+AI  * ZET(2) 

IF  (K.GT.IEND/3)ZETAY=2ET(3)+AI*ZET(4) 

IF  (K.GT.2*IEND/3)  ZETAY=ZET(5)+AI*ZET(6) 


Radial  Modes 

If  the  function  to  be  evaluated  is  not  equal  to  a constant  (one) 
at  the  entrance,  the  pressure  at  those  points  must  be  inputed.  Several 
changes  are  required  in  the  main  program  and  the  SETUP  subroutines  in 
order  to  modify  the  program  for  this  case. 

The  following  lines  should  be  inserted  after  line  28  in  DUCT, 


READ  993  (Y(I),  1=1,  JEND) 

993  FORMAT  (12F6.0) 

PRINT  994  (Y(I),  1=1,  JEND) 

994  FORMAT  (IX,  16F8.3) 
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In  addition,  line  47  in  SETUP2  and  line  51  in  SETUPl  should  be 
replaced  with 

10  Y(I)  = - AA(I)  * Y(I). 

In  SETUP,  line  60  should  be  replaced  by 
10  Y(I)  = - A * Y{I). 


Computations  for  ZETAX  for  a Nonuni  form  Duct 


Recall  that 


’27rnRi 


P * 


and 


p _ -2irniP 
R " Z„ 


where  n is  the  frequency  and  R = ,5e 


-ALPHA -XMAX 


Therefore 


7 -21:711  P 

Pr 


-2irr|ie 


■2irriiR 


- 21771  Rl 


{p  + 2l77li  ) 


(Zuri)^  + 2i7riiR 

1 + {277tiR)^ 


OUTPUT 

The  output  from  DUCT  consists  of  the  following: 

1)  The  solutions  to  the  wave  equation  at  the  indicated  points. 

Note  that  in  reading  the  printout  from  left  to  right,  the  y-coordi- 
nate  varies  while  the  x-coordinate  is  held  constant.  The  last 


20 


column  of  numbers  corresponds  to  the  x-coordinate  of  the  points. 

The  headings  for  the  program  have  been  set  up  for  printing  a 
maximum  of  fifteen  points  in  the  y-di recti  on.  If  desired,  the 
pressure  can  be  printed  at  addition  points,  but  the  headings  will 
be  distorted. 

2)  The  computed  values  for  the  acoustic  power  at  the  entrance  and 
exit  of  the  duct  and  the  resultant  attenuation, 

3)  If  ZXMIN  is  called  and  a multisectional  lining  is  used,  the 

values  of  ZET  will  be  printed  before  the  headings  for  each  iteration. 
When  all  the  iterations  have  been  performed,  the  minimizing  value  of 
ZETAY  for  the  first  lining,  the  attenuation,  the  convergence  criterion, 
the  number  of  iterations,  and  lER,  an  error  condition  will  be 
printed. 

lER  = 0 means  no  errors  occurred  and  convergence  has  been  achieved, 
lER  = 130  means  the  iteration  was  terminated  due  to  excessive  rounding 
errors. 

lER  = 131  means  the  iteration  was  terminated  because  NUMMAX  has  been 
reached. 

Section  VI 
SAMPLE  RUNS 

In  this  section,  the  input  and  output  for  a variety  of  sample 
cases  are  presented. 

Case  I 

Cylindrical  duct:  no  flow,  uniform  lining,  pressure  printed 
at  every  other  point  in  both  directions. 
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Case  I 


Input 


1 


0. 

1. 

.7 

2. 

.5 

50  20 
1 

2 2 20 


Case  II 

Conical  duct:  no  flow,  uniform  lining,  pressure  printed  for 
first  15  points  in  y-di recti on. 

Case  II  Input 

-1 

0 

.9641  .1861 

.3  -.5 

1.  1. 

1 . 

49  20 
1 

2 111 
.75 

Case  III 

Rectangular  duct:  no  flow,  radial  modes,  pressure  printed  for 
first  15  points  in  y-di recti on. 
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Input 


.332 

-.7442 


Case  III 

0 

0. 

1 .0453 
.9287 
1 . 

1 . 

50  20 

1 

2 2 20 
-.775 
.950 
.170 


-1.18  -.440 
-.075  .990 

-.820  -.110 


-1.11  ... 
-.020  ... 
.990  ... 


Case  IV 

Cylindrical  duct:  no  flow,  uniform  lining, 
find  minimizing  value  of  ZETAY. 


Case  IV 

Input 

1 

0. 

1. 

0. 

.6 

-1 .5 

2. 

1 . 

1.25 

60  20 
0 15  2 

0 


Case  V 

Rectangular  Duct:  uniform  flow,  hard  wall 
at  every  other  point  in  y-di recti on. 


Pressures  for  the 
20  nodes  at  the 
entrance 

ZXMIN  called  to 


pressure  printed 
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Case  V Input 

0 

.2 

1.  0. 

999999999999  (20  times) 

1.  1. 

1 . 

25  20 
1 

1 2 20 
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ACOUSTIC  POWER  AT  ENTRANCE  = 1.31163 

ACOUSTIC  POWER  AT  EXIT  = .003C1 

ATTENUATION  = -26,38?35 
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ACOUSTIC  POWER  AT  ENTRANCE  = .65421 

ACOUSTIC  POWER  AT  EXIT  = .02541 

ATTENUATION  = -14.10767 
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ACOUSTIC  POWER  AT  ENTRANCE  = .07034 

ACOUSTIC  POWER  AT  EXIT  = .00193 

ATTENUATION  = -15,  60  759 
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ACOUSTIC  POWER  AT  ENTRANCE  = 2.66222 

ACOUSTIC  POWER  AT  EXIT  = ,64423 

ATTENUATION  = -6.16201 


CYLINDRICAL  DUCT 

0X=  .0328  OY=  ,04?b  ZETAY=  .6000+1-1.5000  ZET4X=  l.OOOO+I  0.0000  £TA=  1.26GI  i=  O.COfl 


o 

C3 

o 

o 

O 

* 

• 

a 

II 

It 

j: 

t: 

Lj 

o 

iS\ 

u% 

{NJ 

C\J 

a 

r-l 

It 

It 

*- 

1— 

Ul 

o 

O 

C5 

O 

O 

o 

i3 

o 

a 

* 

O 

o 

M 

t“l 

o 

o 

o 

o 

a 

o 

a 

a 

a 

II 

II 

X 

X 

aJ 

1- 

1— 

LiJ 

UJ 

fsJ 

CD 

iTk 

o 

CM 

o 

<£> 

U^ 

rsj 

* 

CM 

a 

CM 

(M 

CM 

1 

(M 

t 

l-H 

vD 

M 

po 

4- 

MS 

fO 

-1- 

■ CM 

C» 

a 

CM 

a- 

CM  ^ 

O 

CM 

J- 

CM 

c? 

tH 

u> 

U) 

o 

a 

a 

• 

a 

II 

It 

It 

It 

UJ 

>- 

UJ 

>- 

o 

o 

I t 

T-i 

h“ 

z 

II 

a-l 

1- 

o 

Ul 

UJ 

ti:  i- 

CM 

or 

1— 

CM 

M 

h-  »-4 

M 

lD 

^ X 

z 

X 

aH 

UJ  Hi 

a 

UJ 

UJ 

« 

tO 

iD 

u> 

t— 

1 

1— 

I— 

• 

1—  o 

<c 

H o 

O • 

o • 

c^  a: 

3 

or 

q: 

3 

UJ  UJ 

II 

O II 

UJ 

UJ 

II 

O 11 

-X  :x 

>- 

X 

X 

>" 

o o 

Z 

-I  a 

o 

o 

Z 

-i  CD 

a.  0. 

o 

<£ 

□- 

a 

o 

M 

O 't) 

M 

O cO 

o o 

l-H  CM 

o 

Q 1- 

M CM 

IH  M 

< 

Q£  Kt 

i-l 

M <5 

a:  ro 

h-  1- 

3 

O O 

3 

Q O 

00  (/) 

Z 

z • 

00 

00 

Z 

z • 

3 3 

UJ 

M 

3 

3 

UJ 

O O 

J M 

o 

O 

•- 

-J  II 

U tJ 

>-  X 

o 

O H- 

>*  X 

<t 

o o 

<S 

< 

O a 

32 


ACOUSTIC  POHER  AT  ENTRANCE  = 2.60974 

ACOUSTIC  POHER  AT  EXIT  = .67955 

ATTENUATION  = -5.04377 
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ACOUSTIC  POWER  AT  ENTRANCE  = 2.65959 

ACOUSTIC  POWER  AT  EXIT  = .64267 

ATTENUATION  = -6. 16627 
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ACOUSTIC  POWER  AT  ENTRANCE  = 2.65797 

ACOUSTIC  POWER  AT  EXIT  = .64199 

ATTENUATION  = -6.17:25 


25  20 

DO  rou  WANT  complete  output? 

ENTER  lO'JT,  JOUT,IT/3IJ 


<M 
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I I I I t I 
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K lA  CD  ro  J-  ^ 
CO  hw  ^ 4^  CM  o 


I I I I I I 


CO  CM  <M  CM  O O 
LA  LA  PO  O'  CM  CM  <0 
tH  PO  lA  O)  CO  O'  O' 


PQcdPOO'O'N-'^  *p4POOCOIACMOPOIA  POPO^J"  COCMPOCMCMOO 
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ro  o PO  O'  O'  'O 
PO  PO  oi  vD  LA  CM 
• N1  ^ o CM  vO 


M 1111 


'H  PO  O OO  lA  CM 
® ^ O'  O'  Lf  1 
OO  O'  O'  af>  CA 


III  III 


o po  LA  ro  po  4' 

lA  O PO  4 4* 

^ \D  4 CM  O 


I I I I I t 


C30  CM  PO  CM  CM  O O 
la  lA  PO  O'  CM  CM  ^ 
tH  PO  A M)  CO  O'  O' 


4 oPOO'O'P^^^POooOtACM  aPOlA  POP^4:OCMPO  CMCMOca 
CMPOP^M:>»AArM^cO^£)CJ'0'lAP^lAoPO  44tALAP^CT'CMCM<0 
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096 


IMAGINARY  PART 


xD  ru  c j-  ^ <.0  ^ ^ fn  cj  f\i  o -1-  Ln  liD  rp  Lr>  h "o 

T*  CP  JO  o :o  1^  \0  u>  A t -J"  iO  r^>  ro  CM  t\j  -♦  ^ j 


11 
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ACOUSTIC  POWER  AT  ENTRANCE  = 1.12697 

ACOUSTIC  POWER  AT  EXIT  = 1.14235 

ATTENUATION  = .05309 


GOVERNING  EQ. 

_a!E  + ^+(2irX)*  P*0 

TRANSFORMS  VIA 
•*  •* 

•n  — siny 


TO 


AND  THE  GOVERNN6  EQ. 
BECOMES 


Je 

3X« 


dY* 


+t2irX)*F(X,Y)P»0 


WHERE  FOR  THE 
HYBERBOLC  HORN 

e2X-g-gx^2(cos^y-sln^y) 
F{X,Y)»  - 


4Y 


Example  of  conformal  map. 


FIGURE  1 


37 


Sound 

Source 


Conical  Duct 


Sound 

Source 


Catenoidal  Duct 


Sound 

Source 


Examples  of  Variable  Area  Ducts 
FIGURE  2 
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FIGURE  3 
39 


APPENDIX 


program  duck  input, output,  TAPf  1,  TAPE:  2,  TAPt  3,  TaPc  4 ) 
COMMON/3 IG/AI ,2ETAX,ZETAr ,XMAX,yMAX,£TA ,XM, I £ ND , JEN C , uY JX , M , NN , 
ICON,  A,B33e,a(20,20)  ,C(  20,2  J)  ,T  (2C,2C)  ,O,Y(20)  ,IP(2:)  , V(2j)  ,JX,DY, 
2 AAt2(iKCI  IT  (20  » , lOUT  , JOJT,  KOUT,  IT  , IDUCT  , KEN  J ,Dd 
COMPLEX  3,C ,T, 0, Y,V, AI,ZETAX,Z£TAY  ,C1,  A,Al,bl,AA 
external  FUNCT 

DIMENSION  ZET ( 2) , G(2) , HM( , NA (6) 

*1  = 20 


931 


PRINT  961 
read  963,  lOUCT 
FORMAT(*  enter  DUCT  GEOM\ 


PRINT  991 
read  990, 
PRINT  983 
READ  990, 
PRINT  987 
READ  990, 
PRINT  968 


XM 

ZETAX 

ZETAY 


NO MUNI FORM, REG T ,CYlIND;-1,l,1 


READ  990,  XMAX,  YMAX 
PRINT  985 
READ  990,  ETA 
PRINT  981. 

READ  992, lEND, JEND 
PRINT  992, lEND, JENO 
PRINT  962 

READ  992,KOUT,NUMMAX,IH 
PRINT  969 

READ  992,  IOUT,JOUT,lT 
989  FORMAT (IX ,*£NTER  lOUT , JOUT , I T/31i* ) 
992  FORMAKZIS) 


982  FORMAK*  do  YOU  WANT  COMPLETE  OUTPUT?  T OR  1 *) 

983  FORMAT (313, F10,0) 

984  FORMAK*  enter  I END  ( X ) , J END  ( YJ  213 

985  FORMATt*  ENTER  ETA  F10.3  ♦) 

98b  FORMAK*  ENTER  XMAX, YMAX  F IC  , i,  ,F  1 C . C 

987  FORMAK*  ENTER  ZETAY  FlO.O.FlO.u  ♦) 

988  FORMAK*  ENTER  ZETAX  F13.0,F10.0  *) 

990  FORMATCaFlO .0) 

991  FORMAT (*  ENTER  MACH  NO.  FIO.O 
A1=(D.,1*) 

ZETAX=ZETAX^(1 .*XHJ 


NN=JEND 


DX  = XMAX/(IEND+  1) 

OV  = YHAX/(  JENO+1) 

DY DX=OY/aX 

C0N=2.»3. 14159265358979*ETA 
1 A =-(  l.-XM»XM)  *DY0X*»2-C0N*XM*AI*DYL>X*0y 
0a8a=2.+2.*DYDX»*2-CON**2*DY**2 
BB39=0B8B-2 .*X XH* D YD X** 2 
ZET  (1)=REAL(ZETAY) 

ZET  (2)=AIMAG(ZETAY) 


40 


IF (NUMMAX . GT. 3)  GO  TO  ICO 
GALL  FUNGT (2,ZeT,FMIN) 

GO  TO  f03 
IC:  NVAR-2 
IOPT=0 

CALL  ZXHlN<FUi4GT,NVAR,  H , tJUM^lAX,  I OPT  , ZFT  ,HH,  G , FMI  N , h A , 1 Ek) 

PRINT  gyg, (ZET (I) ,I=1»2> (F^IN, IH,NUMMAX, lER 
q99  FORMAT (IX ,3F16 .8,315) 

50C  STOP 
ENO 
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o o 


sufi^nijTiN'^  S'-'^up 


r 

n C'^LTMPPTr:4L  NON-UNI^’OP^I  HUQT  USING  CONPOP'^fiL 


PQHMON/  UG/  ai,  ZLTAX,  ^e-T-'r  ,yMtX,  v^A  X,£T  a , X'^  ,IlNO  t J- NO  t D VOX  , M , NN  t 
iPGM , a , ogfji,  H » t c ( 2 ,t (2„  ,2  , j »n, v (2  »,  ip  ( 2w ),  v(  nx.ovt 

2 A4  (2.:)  ,GUT(  2:»  , IP'JT,  JOUT,KDUT,IT,nuOT,KPNn,OQ 
COHPLPX  P.CtT.TtY  ,V,  A-»7ETflX,  ZcTAY  ,C1»  .A,ai,31,Afl 
PPWIMP  1 
P^INr  ggo 

peat  alpha 

qQ%  FOPNAT (2P1„ . 0 > 

ggg  F0P''AT(1X,*  E.MT£P  ALPHA, FlO.^  *► 

A=-ornx*»?*aLPHa»nY*nx'Dxx?, 

PH  101  J=lt  I'HD 
i':r  Aa(  n = A 

00  f=.  '<=1,TP*I0 
PO  T=l«.t-Nn 

P ( I , I > = -OYP  <*•’- ALPHA*nY*nYDy/2. 

OP  5 J=i,JEHn 
TF(J.U‘.II  0(  T,J>=1. 

P ( I ♦ 1)  = A . 

YK=  , 5 

Apr,l  = OX*FLOAT  <KI 
AR02=aLPHA*FLPAT( T)*0Y 
FOPY=-nOS(  ARG21 / SIN ( SPG?) 

SHA'srsALPHA*  ALPHA  ♦ExP  (2.  •alpha  *AFGl)  •XK*XK 
TF(I.^0.J»  ■>  ( r,  J J =d330  ♦ { 1,-SHAPEI  *caN**?*nY*»2 
IPfJ  .EO.  T-1)  HI  T,  J)  =-l.-ALPHa*F0PY*0V/2. 

IP(  I .EO.  I*-1T  0(  I , i»  =- 1.  • ALPHa»F0PV*Cx/2. 

^ f'ONTTNxr 


P TNCLUP"  lOUNOA-v  gATA  POEPESPONPTNG  TO  Y = '‘  AND  Y = 1 

Pl  = -?. •CON*AI*nY*  ALPHn»yK*EXP (ALPHA^ ARGl ) XZ'TAY 
B(  JENO, JEMn-ll  = -?. 

P(l.  21  = -'’  . 

AKG2=ALPHa»'-LPAT  I J"NO>  •py 
FnFV=-C0S(APG2»/SIN(APG^t 

Pt  JENP,  J-NO)  =P(  JEHQ,  J-NOU- (ALpHA*F0FY*0Y/2,  - l.J'Rl 

IF  ( K .EO.  T - NP)  GO  TO  G 
WPITEdl  AA.n.O 
0 rONTINU'^ 

TNOLUOe  SOUNOA'Y  OATA  C Q^P ESP CND TNG  TO  X=1 


shapf=alpha»xk*exp( alpha^argi ) 

CQNl=:nON 

CON=CON»SHAPF 

DO  7 T=l,JENn 

7 0(I,I)  = B(I.:t-C(I.I)  *7. •OON*AI*aX/ZETAX 
CON=CONl 
WPITFdt  AA,B.r 
CO  <1  J = l,JPNO 
3 AA(  l>=A 

DO  1C  T=1.J"ND 
10  Ytrt=-A 

PFWTNO  1 .p 

RETUPN 

EMC 


SUBROUTINE  BLOCK 

CO*1HON/3IG/AT,ZETAX*7ETflY,XMAX,YHfiX,ETfl,X‘4,lEND,  JEND*DTOX,M,NNt 
ICON,  A,BBRBt3(20  f 20  »C  (ZDiZOI  , T (2C * 2Q  ) ,D , Y (2C’J  tlP  ( 2«)  , V t 23 ) , 0X,0 Y, 
2 AA ! 20» ,OUT<  2C» , I OUT , JOUT , KOUT ,IT , I DUCT .K END, 08 
COMOLEX  B,C»T ,D, Y,V, AI ♦ ZET AX, ZFTAY  tCl,  A , A1 , Bl , AA , CC 
OTHENSION  CCI2D,2r) 

PTMENSION  TN0EX(3*J1» 

CALL  0PENMSC2, INDEX, 311,0 
REMIND  1 
P^MINO  3 
PFWTNO  u 

PEA0(4>  (V(I», l=l,JENn) 

REAOtH  AA,D,CC 

call  NRITHS<2,n,8r-‘’,l,-l,'') 

CALL  WPITMS(2,CC, 800 ,2,-1, D) 

CALL  MRITHSf ?,Y ,4r, 

OP  100  K = ’,'^ENO 
DO  12  1 = 1, .rND 

12  V<TI=Y(Tj 

REao{4)  <Y( 1=1, JEND> 

DO  1*=^  T = 1,J-M0 
no  15  J=l,J-Nn 
c(i, J>=cn (I, J> 

15  T(I , J)=B(T,JI 

call  LU3(NN,T ,M,TPJ 
CALL  L='OSt<NN,T,m»IP,  V) 
no  2r  i=i,j'^NO 
A=Aa  (H 

TF(K  .EO.  T'Nni  fl  = A + C(T,T) 

D=v  (T) 

21  v<it=o-a»v(Tt 
pEAodi  aA,D,nc 
no  50  T=1,JENP 
np  T-d  J=1,JEMD 
■'ll  V(J1=C(J,T1 

call  LEQS3<NM,T  V) 

pp  UK'  J=l,J£Nn 
A=AA  t J1 

TF{K  .EO.  I'ND)  A = A + C<J,J> 

40  Dt  J,T»=Bt  J.n -A*V  ( J1 
= CONTTNU'' 

kk=  >*k-t 

CALL  WOIT'^St  2 ,P,8G  1,  Kif  + lt-lir  ) 

CALL  wpiTHS (?,r, pro, <*^*2,-1,^ ) 
call  WpITms  f 2,v,4C  ,i<K«-3,-l,n 
IP'’  CONTTNU'^ 
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O 


INVENT  T^'IflNGULAR  MATPIX  TO  OBTAIN  SOLUTION 

no  iin  T=i.JFNn 

V(I) =V(T» 

nn  in  j=i,jeno 

113  Tf I , J»=R(T, J| 

TALL  LU3 ( NN, T ,M,TP> 

TALL  L=’OST(NN,T,M,IP,  V) 
no  IHu  I = l,.>FNO 
12^  Y(TV=V(T) 

P‘‘WTND  1 

WPITe('>  (Y ( I ) , 1= 1 , J€ND» 

TrM31-=T!^Nn-l 

no  K^ltlENOl 

KK=  j»T- N''-"*K-'' 

CALL  9E AHNSt  2,B,  ».%K'<'«-1) 

CALL  ( -»  ,C  , «<■  0,  KXt  2) 

call  ^"A0MS(  2,AA,4s.,  K<«-T> 

nn  1 4u  T = i^jrN|p 

V(  . 

n*'  1 ■'J  1=  1,  ) run 

1',  vfIt=CtT,J)*v(j)+v<I) 

ILn  Vf  f ) = AA  ( 1 1-v/ ( I> 

PALL  LU-ff  NN,'3,M,I  P) 
nALL  L- nS*'(  NN,D  ,N  ,T», 'M 

no  IfiL  T = l,  )PNn 
lei  Y(  :»='/(  T) 

WFT'rr(^)  (V  ( I ) ,T  = 1 ,.n‘'n) 

2r''  rON'^i'j'JE 
c-WTfgn  1 
CALL  CLOnM^'CEt 
o"T'loN 
'^N'l 
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SUBROUTINE  OUTPUT 

COMMON/ 3 IS/ A I jZETAX,  ZETAY,  XMAX,YMAX  ,ETA  ,Xi1,  I£t\lO  , JtND  , DYDX  ,1,  NN, 
10  ONtAtSOBa,  6(20,2  0)  , C ( 20  1 2 Q>  » T ( 2C  , 20  ) i D i Y <20 ) , IP  (2L- ) , V (2  j ) , JX  ,ljY  , 
2 A A (20), OUT (20 ) , lOUT , JOUT, ^OUT, IT , IDUCI , nEND,03 
complex  e,C ,T,0,Y,tf , AI  ,ZETAX,ZETAY  ,C1,  A,Al,ai,AA 
PRINT  996 

DATA  YEQ/3H  Y=/,XEQ/3H  X=/ 

IFdOUCT  .GT.  3)  PRINT  993 
IFdDUCT  .EG.  0 ) PRINT  9^2 
IFdOUCT  .LT.  0 ) PRINT  991 
D91  FORMAT  (*  NON-UNIFORM  DUCT  '^) 

992  FORMAK*  RECTANGULAR  DUCT  *) 

993  FORMAK*  CYLINDRICAL  DUCT  *) 

PRINT  995,  DX,  OY,  ZETAY,  ZETAX,  tTA  , XM 
IFCKOUT  .Ea.  3)  RETURN 
PRINT  989 
REWIND  3 

DO  200  J=l, IT, JOUT 
20C  OUT(J)=J*DY 

PRINT  990, ( (YEQ,OUT (J) ) , J=1,IT , JOUT) , XEG 
990  FORMAT  (IX,  15  (IX,  A3,  F-i.  2)  ,3X,A3) 

L=0 

DO  210  I=1,IEND 
REAO(3)  (Y ( J) , J=l, JEND) 

L=L  + 1 

IF(L  -LT.  lour  .AND.  I .GT.  1)  GO  TO  210 
DO  209  J=l, JEND, JOUT 

209  OUT(  J)=REAL  (Y{  J)  ) 

X= (IENO-I*l)*DX 

PRINT  999, (OUT( J) , J=1,IT, JOUT),X 
L = 0 

210  CONTINUE 
REWIND  3 
PRINT  999 
PRINT  995 
PRINT  998 

DO  201  J=l,IT,JOUT 
201  OUT(J)=J»DY 

PRINT  99  0, ( (YEa,OUT ( J) ) , J= 1, IT , JOUT) , XEQ 
L=0 

DO  212  I=1,1EN0 
REA0(3)  (Y ( J) , J=1,JEND) 

L=Lfrl 

IF(L  .LT.  lOUT  .AND.  I 
DO  211  J=1 , JENO, JOUT 

211  OUTt J>=AIMAG(Y  (J) ) 

X=  dEND-I  + l)*OX 

PRINT  999, (OUT(J»,J=l,IT, JOUT) ,X 
L=0 

212  CONTINUE 


,GT.  1)  GO  TO  212 


45 


994  FORMAT (IX, l2Fia.5) 

995  FORMAK  23X,  ♦IMAGINARY  PART^) 

FORMAT  (IX,  ♦LX=»,F6.4,»  OY  = »,Fb.H,*  Z ET  A Y = * , F 7 . 4 , ♦ ♦ I • , F 7 
1 • ZETAX  = * ,F7  .4,»  + I»  , F7.‘-,»  ETa  = *,F7.4,*  m = *,F7.4,//) 
989  FORMAT (21X , ♦REAL  PART^,//) 

997  FORMAT(Hl) 

998  FORMAT (IX, 7//) 

P99  F0RMAT(1X,1EF8.3) 

RETURN 

END 


• » 
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rutl  ROUTT  M?:  ';OUNO 

rOMMON/OIGi^AI,  ZETflX  , 7E T AY  , XM  AX,  YMfl  Y , ET  A , XM  , ItNO  , JHNO  ,0  YDX  , H,  NN  , 
ICON,  A,BR0'?,  0(2C,2CI  ,C(20,20»  ,TI2C  ,2CI  ,n,Y  (20y  J ,Vt?'5)  ,DX,DY, 

2 AA (231 ,0UT(2P» , TOUT, JOUT , KOUT , IT , lOUCT ,KFNn, OB 
COMPLEX  R,C,T,0,Y,V,AI, ZETAX,ZETAY  ,C1,  A,A1,31,AA 

PfMTNn  7 
■DA  = nY 

IF  ( XM  .ME.  0.  ) GO  TO  299 
SUM  1=0. 

PFAn(T»  tYUI  ,T=1, JENO) 

(Vdt  ,1  = 1 .J'^NO) 

no  220  J=l,JPNn 

IFtTOUCT  .NE.  P)  nA=J»OY 
220  SUM1=SUM1+REAH  AI*C0NJG(T(  J)»  *(Y(  J|  -VtJl  t » 

1 

SUM2=l, 

KLaST=TEND-7 
DO  225  K=l,KLa^T 
REa0t3)  (Yd)  ,1=1  , JENO) 

2^S  rONTTMUE 

READd)  (Vtl)  ,T  = l,JENO) 
r'O  230  J = l,  JFMH 
TF(TOUCT  .N",  D na=j*0Y 

<;U'^Z  = SI)M2fR- AL  ( AI^CONJG  tV(J>)*(Y(J)-Y(J)y) 

1 •Ofl 

PPINT  99« 

PPTNT  9P4  ,SUM2,SU'‘'l 

IF«  sum2.lt.  ;. ) c'^turn 

no=l'.*ALOG10 (SUM 1/SUM 2 » 

PRINT  790,00 
Ri^TURM 

297  I^N72  = Tp.Nn-'' 

Pi^AOtSl  <V(  J t ,J=1 , JFN'') 

PFA0(3»  (V(  JJ  ,J=1  , 

x=  XMflX-OX 
no  300  J=1,J'^MD 
T ( 1 , J » = T ( J ) 

3C ''  Bt  1 , jy  =cExn  (coN^a  i^x/vM)  •(Y(j)-v(jyy/r)x 
no  aCO  T=l,TEMn2 
FacT=FLnaT( I ) -2 , * float ( i/zy 
WT  = 4. 

IF(F1CT  .LT.  .sy  NT=2, 
y=  XM AX-OX- TX^FLOATt  T) 

Pt'A0(31  < AA  ( jy  , J=  1,  JENOy 

DO  ^5*’  J-1,  l-NH 

ci= ( Yt jy -AA  t ji y / ( 2. ■OX) 

n(l,jy=9»l,J>*WT*nt*c~X3(CCM'‘nT*-X/XM  ) 

Y ( j I = V ( j y 
■''^y  v(jy  = AA(jy 
CONTTNU'^ 
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nr  4 5C>  j = i,j^ND 
r*!  = ( V ( .1 ) - fl  fi  ( j n / n y 

T(i, j»=n(i,  j>  *c£vo(CnN*fti*nx/yM) *01 
V ( J » = ^ r 1 ♦ J J 
ni  = ’^x*g(  1,  j»  /?, 

x=  Ax-ny 

ri=noN*Y(j)*c"yo(f'0’'i*Ai*x/yMi 
n(i , j » = ( AtfCii*CE  y’c  - ''0N*Ai*ny/yM» 
t.5''  rONTTM'J^ 

ri=  CON*  A i/(  y '(•xM't-vM) 

X=yMay_gY 

ni3-xM*r^xD(ci*x)+XM*''ryp(ci*nx) 

IP  ( WT.PO.  ) ppImt  qq-r 

qo3  roq‘iaT(.  liTMO  ir  rycKj  yTELOING  TN  nUfiDKA TU^E* ) 

"p  Ri*'  j = i,  iPMn 

TP(l,c'?.l»  n = (’^(  J*i  » -Y(  J»  I /pY 
Tr(  J .CO.  J-NU  A 1 = ( V ( J»-Y  (J-1)  )/'’Y 

’PC)  .GT.  1 .and*  j .lt.  jen^)  ai= { y ( j*i ) -y { j-1 I » / ( ?, *py > 

Yi  = fi.fy'i)*noNjr, (Y(j))*Y(j>  xM*(i.«.xMt*(i. fXN) *ai*coNjr,( fli)/ 

1 ( . *r;'^N*r<^«n 

IP(  Y-.unT  .^’P.  D nA  = J*OY 
SUM  = +''  “I  * XT 

rp''  ^ON^’INUp 
p|)M  1 = '’  . 

HP  ftp"  j=i,>jpwn 

YF(J.C').:)  51  = (Afl(J«-l)-AA(jn  /UY 
TF(J  .CO.  JCNC>  A { AA  ( J)  - AA  ( J-1 ) ) /OY 

IP(  I .';T.  1 .AMO.  J ,CT.  JENO)  Al=(AA(J*l)-Aft(J-in  /(E.*nv) 

yptj  .pp.  11  E!l=  ( P(l,  ’J -fU  1,1  M/OY 

TPtvJ  .CO.  JCN^t  01=  (3 ' 1,  JPNO)  -Rd  ,JtNO- 1)  )/nY 

TP(J  .Pd.  1 .am 3.  J .lt.  jeNH)  0i  = (T(i,J«-1>-3(1,J-H)/1’.*OY1 

Cl  = AI*Al  + AT*  X'^*  31 /CON 

yi  = ^EAL  (CO’^'JO  (AA<  JM  (l,.l)/rrN> 

XT  = yT*xy*''Of;jr.<AA(jn*AA(j)/?. 

yT  = YT4-yM*  (C'^N  jr,  (q  ( 1,  J ) ) ( 1,  J ) ^-coNjr.f  oil  * ri ) / (2.*C0N**C) 

IPf  TOUCT  .\'P.  "1  OA  = J*nv 

SUN1  =S'JN1  +01*  XT 
f.  • ■:  r'^fdlNM' 

DtriMT  qqn 

PPINT  334  , SiJMi,<:uM2 
IPfSU’n  .LY.  .)  rs'Tu^N 
r»p=  1ft , * ALOG 1 p (SU'^^/su‘'i  y 

C^I'JT  ^QT.ori 

aqt,  FnPMaTi  i y ,* ->  coLdTTO  OPNCP  AT  rNTPANTC  = *,PL  .5,/. 

1 1 y . * AC'iUSTTC  P0WF3  A7  CyTT  = *,p:'!.pi 
9^17  pn?  way  ( IH  1) 

OG3  po^NAT ( 1 Y, ///) 

gnq  FP?M  A t ( 1 y , y A T T^  MU  A TT  0 N = *,pi'.3) 

P^TUPN 

CM  q 
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c_'  c l:  c cj  c: 


SUnoO'JTTNt 


TMIS  'JUG'^niJTTMf  S"T<:-Uo  THT  PINIT£  OlFP-£Ri£Nnc  ’^air^TX  PQK  A -’cTTANbULAR 

•“O^^MnM/lTG/AI  » 7f  TAX  , T4Y  tX^*  tX,v'iAX,  "Tfl  , X^^IEMn  , JE  ><0  « 0 '^HX  , M , NN  , 
tnOM^A.GBa^**  1(2f  ,2'*  > ,Ct2.:,  23)  » T(a.  ,2«»  Y(3'  j ,IP<^.:>  ,V(E:)  ,OX,Oy» 

2 4A(2'')  ,nUT(2'’)  ,TOU’’»iOUT,<OUT,IT  , I OUHT  * K £N  Hf  08 
CO’-IOLEX  '^tCtT  ,G»  Y ,V,  AT  , 2ET6X,  7rTAY  tOlt  A,A1,01,AA 
PEMTMn  1 
DO  6 K=l.TENin 
r>0  5 T = l»Jc‘'D 

C(  I,T>  rfi*-?,  *CON*XM*AI  •'’YCX*nY 

nn  c J=l,.ieM0 
irij.ME.T)  D(  I,J)=3. 
n ( I , J ) = " * 

IF(I.EQ.J)  K T,J)  = 3'‘98 
TF ( J ,E0. I+ir  B(  I,J ) =-l. 

TFtJ.El.I-1)  B(  ItJ)=-l. 

B CONTINUE 

IMrLUT=’  “DlJ'IOARY  OATfl  r 0?R  =SP ONDTNG  TO  Y = C AND  y = i 
ni=  7,  *X’^*X'^*AT*nY0X/(C0’^*7FT  AY*nX) 

Pl=-2«*  A1-2.  *C0N’ AI’Oxy^ETAV 
C1 = Al-2«*XH*nY0X/7-T AY 
Al=  Al*'> . •XH'‘nYnX/7£T  AY 
R(  1,2)=-?. 

R(  .lEND,  JEflD-1)  =-?. 

B(i,n  = Gn*  i»-ni 

R(  JEND,  =B(  JENDt  JEND)  - 81 

r(i,i) =C(t, 1) -ri 
0(  JENO,  J'^Nn)  =C(  JENDI  - 01 

TFtK  .ED.  TENDI  GO  TO  8 
Pf'  10  3 J=1,JEWD 
iCr  AA(  )»=A 

AAf II =AA  1 1) - A1 
AA  ( lENT)  = AA  ( JEND)  -A1 
W=»ITEHl  AA.Q.n 
f,  nONTINI)'^ 
n 

C INOLUO-  BOUNHARY  DATA  GORR'^*i”ONni NG  TO  X = 1 

r 

no  7 l=i,JENO 

7 R( I. T)  = B(T, II -0 ( r ,11  *7,*00N*A  7*DX/Z- TAX 
PO  1 J=l,  0*^8  0 

8 AAtl)  = A 

AAt  1)  =AA  ( 1)  -A1 
AA  (J<^NO)  =AA  ( J«"NDI  -A1 
WPITCdl  AA,B,C 
DO  10  I=1,J'N0 
1"  YtI1=-AA(I) 

PFWINO  1 
RETURN 

FMH 
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O c c Cl 


<;U‘1R0UTTNE  SETUP2 


n 

THIS  SUBROUTINE  SETS-UP  THE  FINITE  DIFFERENCE  MaTRIX  FOR  A CVLINDkICAL 


rONMQN/RIG/AT, ZFTAX ,7ETaY ,XMAXtYHAX,FTa, XM«IENDt J-ND,D¥DX, M, NN » 
irON,a,q03B,3(2O»2C»,C(2Uf  2C1 tT(2C  t2At ,n,Y(2„l  tIPt23) ,V(23> ,DX,DY. 
2 AA (201 ,0UT{2"»  » I OUT, JOUT , KOUT , IT , I DUCT » KEND , OB 
COMPLEX  i3,C,T  ,C,  Y , AI  , Z-TAX,  Z«^TAY  ,Cl,  a,Al,ai,Aa 
R'^WINO  1 
DO  6 K=l,TFMn 
DO  5 1=1, JEN n 

Cf  1,11 = a+2 . *CON*XN*ai*nYDX*DV 


DO  5 J=1,JEND 

TFf  J.NE.T)  0(  = 

R(  I,J1=P. 

TFd.ED.Jl  n(  I,J)  = B"fl'^ 
TF{J  .EQ.  I+l) 

IFtJ  .EQ.  T-1) 


B(I,Jt=-l.-l. /FLOAT (2*1) 
B(I,Jl=-i.+l./FLOAT(Z*Il 


S  rONTINUE 


TNCLUD-  nOUNDA^v'  naiA  CO'^’PESDONOING  TO  Y='u  AND  v = l 

ai=  2.  '^XM»XM*AT*"tYQX  / (CnN^ZETAY^nXl 
Pl=-’.*Al-2. •COM*ai* QV/ZETAV 
,*XM*DYDX/ZFTaY 
ai=  AH-2  .•XM*ovnx/Z£T  AV 
R(  l,Z)=-2. 

Bf  J'^ND,  JEND-1)  =-2. 

or JENO, JENDl  = B( JFMD, JEND) - /FLO  AT (2* JEND) > 

OCJEND, JEMDl =r( JENO, JENOl - ri*(l.+l,  /FLO AT ( 2* J END) 1 
IFtK  .EQ.  I"ND)  GO  TO  6 
no  inn  j=i,JENT' 

IC''  Aa(J)=A 

AA ( JEND)  = AA ( JENDl - Al • < 1 . +1 . /F LOAT ( 2* J EM  01  1 
WPITtdl  AA.B.C 

6 OONTINUP 
C 

C INCLUDE  ROUN'^ARY  DATA  C ORRE^TP  OND I NG  TO  X = 1 

o 

^0  7 1 = 1,  JEMO 

7 R(I,I)  = '3tI,I)-P(I,I)  *7  . •rf>fj4AT»DX/7ETax 
^0  J = l,l“Nn 

8 aa(ji=A 

aa<.IENQl  = AA{J£NP)-ai*(l.tl.  /float  (2*  J>-MD)  y 
W^ITEd)  AA.S.C 
DO  in  I=l,J£Nn 
ir  v(T)=-aa(ii 
RFWINn  I 
pETijon 
FMD 
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SUn'?OIJTlNE  '^UMPTC  N,ZE:T, 

C0MM0N/3IG/flI, ZFTAX,ZETftY,XMAX,YM4X,ETAtXM,I£NO,JENO,DYOK,H,NN, 
lC0N,A,RR33,3C2t  .2GI  »C(23,2i]t  ,T(2ri,21l  ,0,  Y (20  , IP  ( 20  ) , V(  20  ) , OX.DY, 
AA  ( 20  tOUT(  20  , lOUT,  JOUT,  KOUT  , IT  » T OUCT  , KENn,  OB 
rOMPLEV  B,C,T,D,Y,tf ,flT»7ETAV,ZETAX,Cl, A, A1,B1,AA 
pTM-NSTON  ZET(Mt 
ZFTAY  = 7ET(1)  ♦^AI*ZF7(2) 

IF  (7ET(11  .L'T.C..  J GO  TO  4'' 

A=-  tl,-X*1*XM)  •OYnv*»?-CON»XM*  AT*OYOX*DY 
TF  ( TOU''T,LT.  C I call  S-TIJP 
TnuCT.  £0,  P)  CALL  SETU'^l 
IP ( tOUCT.GT. n ) CALL  SETU^^ 
ot^'JTNO  4 

WRT’'P<4)  (Y(Tt,I=l,  JE'^'O) 

OP  I'-  T = i,J»-Nn 
li]  Y(T)=-'. 

TEN  U-TEMP-1 

no  IE  J=l,I'Nni 

W=‘ITFf41  (V  ( I ) , 1=  JPMn) 

IE  CONTIN'.JE 

WDIt£(41  aa 

CALL  BLO-'K 
call  O'JTPUT 

op=  ■ . 

CALL  SHUN" 

GO 

4"  np=T, 

EO  P=OB 

RPT’RN 

rwn 
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w»«  ««  «««**  *«»*«««  ««»«  **************■¥  «»»*«««  *.»»«»««**«« 

C TOFNTTFinflTTON 

r,  LU'^  - LU  P'flCTnRT7fl'^ION  OF  fl  *?EaL  SQUARE  MATRIX 

c pr*!>Ton(^  ?l>5poutine  subprogram 

r Ae-ROSOACE  ^FFEARCH  LABORATORIES 

C WRIGHT -D A TTEP SOM  ONTO  ‘45413 

C cu=>P0SF 

0 LU3  co“”UTFs  Triangular  matrices  l and  u and  a '»-rmutation 

C matrix  R RATTSFYING  LU  = PA  GT rfEN  AN  N-OOUApE  SEAL  MATRIX  a. 

C LU’  IS  TMTENDFF  fO”'  use  with  THT  entry  point  LTOSJ  tj  produce 

0 SOLUTIONS  OF  THtr  V'CTOF  EQUATION  AX  = R. 

0 rONTOOL 


r DIMENSION  A(M,M),  tp(M) 

r 


C 

r, 

C 

C 

0 

0 

c 

c 

0 

c 

c 

c 

r 

c 

c 

0 


« 

CALL  LUT (N tA, M, tO) 

WH'RP 

N IS  AN  INTEGER  INPUT  VARIABLE,  THE  OROER  0*^  THE  matriy  A. 

A AS  A real  input  ARRAY  IS  THE  MATRTX  TD  SE  TPI  ANGUL  API7tb. 

A AS  A REAL  OUTPU'^  ARRAY  IS  THE  UPPER  TRIANGULAR  FACIOP  U 
TN  0(4-,  jj,  I ,L-",  J,  ANP  THii  LOW="R  TRIANGULAR  FACTOR  t IN 
A 1 1, J)  , T .FT.  J. 

M IF  AN  INT'^GFR  INPUT  VAFIAQLE,  THE  ROW  DIMfNSION  0^  A. 

IP  IS  AN  TNT'^GER  OUTPUT  APFAY,  IP(K),  X .LT.  N,  BEING  THE  IN- 
DEX OF  THE  K-TH  PIVOT  BQW  WITH  IP(N)  BEING  Z“RO  IF  A 15 
SINGULAR  AND  OET(P)  IF  A IS  NONSTNGLLAP, 
other  PROGRAMMING  INRORMATIOM 

LU3  contains  the  FORTRAN  STATEMENT 


C data  NIIM/^cr/ 

C 

C NUN  IS  initial  argument  SUPPLIED  TO  Thf  SYSTEM  SUBROUTINE  cN- 

F TEDED  IN  CASE  LU3  WAS  PHTERED  WITH  A NONPDSITIVE  N. 

F OETtAt  may  be  CALCULATdO  FROM  THE  FORMULA 

C DET  a = FLOAT  tlP(M)  ) • A ( 1,1  J » ...  •A(H,N» 

C LEOSS  IS  CALLED  TO  COMPUTE  SOLUTIONS  OF  LINEAR  SYSTEMS  AFTtR 

c tpiangularizatton  of  a by  LUT. 

0 other  PROGpAMS  required 

c the  system  subroutine  is  called  for  FPROR  tracing  ANO  TERM- 

C INATION. 

r method 

0 THR  MATRIX  A IS  REOUCEO  IN  SITU  TO  TRIANGULAR  PRODUCT  FORM 

C USING  GAUSSIAN  ELIMINATION  WITH  PARTIAL  PIVOTING. 
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O f_  C'  c t c e CJ  C'  C c;  CJ  C L c O O C.'  c o c c c c o o O O CJ  u t:  t C'  O C 


C nENTT^'T^ATTOM 

C LEO"'?  • SOLUTION  OF  4 LINEAR  SYSTEM  GIVE''’  A TRIANGULAR  FACTO® 

r.  TZATIDN  OF  THE  OOEFFiriENT  MATRIX 

r.  FORTRAN  S'Jo='CI!'^TME  SUBPROGRAM 

C AER0S®ACR  RFS£fiDCH  LABORATORIES 

0 WRIGHT-oA  TTFR*^ON  AFR,  OHIO  45t.3'» 

C ou9Pf^<^r 

n LRQS''  C0‘*‘>UT''S  the  SOLUTION  X OF  THF  LINEAR  SYSTEM  LUX  = PAX 

r = pq  HMF-»e  L,  U,  and  ” APE  COHOUTtO  FROM  A QY  LUl . 

'CONTROL 

HIM- NS  ION  A(m,m),  IO(NTt  BIN) 


CALL  L'^0S’(H,A  ,M, 


R,9) 


mhe®'^ 

N IS  AN  TNTr^cp  input  VARTABL 
A IS  A vEAL  input  array,  THE 
COMPUTED  BY  Lin. 

M IS  an  INTEGFP  input  VAFIABL 
IP  IS  AN  IMTEGfR  input  ARRAY, 
COMPiJTED  OY  lU3. 

P AS  A RfAL  IN=UT  ARRAY  IS  TH 
R AS  A f-'EAL  OUTPUT  APRAV  IS  T 
OTH^p  programming  INFORMATION 
LEQS’  TS  an  entry  point  to  SUB 
other  procpahs  RFOUIREO 

NON" 

METHOO 

TH*^  THo  TRIANGULAR  SYSTEMS  LY 
TURN  FOP  THE  SOLUTION  X OF  LUX 
NONSINGULARITY  np  THE  MATRIX  U 
ROGA^ING  IP(N»,  AH  ATTEMPT  TO 
WILL  RE'^ULT  TM  AN  TNFTNITF  CR 


RFFERt^UCES 

tt)  CL‘‘VE  B.  MPLER.  ALGORITHM 
COMM,  ACM  15ll9''2t,  P.274, 
(2)  PAUL  J,  NIXOLAT,  the  ARL  L 
TECHNICAL  FEPORT  Apl  71-iJl 


E*  THE  OFDE"^  OF  THE  MATRIX  A. 
TRIANGUL  ARIZED  COF.FFICTENT  MATRIX 

E«  THE  ROW  DIMENSION  OF  A. 

THE  indexes  OF  THE  PIVOT  ROWS 

E COLUMN  VECTOR  RIGHT-HAND  SIDE  P 
HE  COLUMN  VECTOR  SOLUTION  X. 

routine  LUl. 


= PB  ANO  UX  = Y ARE  SOLVED  IN 
= PB.  NO  CHECK  IS  MADE  FOR 
. THIS  CHECK  IS  MAOE  BY  INTi_R- 
SOLVF  UX  = V WHERE  A TS  SINGULAR 
INDEFINITE  QUOTIENT, 

423  - LINEAR  EQUATION  SOLVER, 

INEAR  ALGEBRA  LIBRARY,  ARL 
T7 (1971)  . 


SUBROUTIN'"  LU3CN,  A,  M,  TP,  B) 


triahGULARIZAtioN  of  a R'"AL  SOUAPE  matrix  using  GAUSSIAN  ELIMINATION 

DIM’^NSION  A(H,M),  IPtN),  B(N) 

COMPLEX  A,  B,  T 
EQUIVALENCF  (KPt,  KMl) 

OATA  NIJM  /23<’/ 

TP(N)  = 1 
NHt  = N - 1 
IF  (NMD  2,  65,  4 
R CALL  SYSTEMtNUM,  12H1ILLFGAL  ARG) 
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o o 


4 on  K = 1,  NMl 
KOI  r K 1 
L = 

00  13  I = KPlt  N 

IF  ( f ABSIRFAH  ft  ( I »K)  H ♦ flBS(AIMftGf  ftdtKI)  n .GT, 

* ( ftBSCRFftL(  A (L  tKn  ) ♦ ABS  ( A I MAG  ( A (L,  K » > M ) L = I 

10  CONTINUF 

IPCO  = L 

IP(N)  = ISIGMd,  K - L)»IP(Nl 
T = AdtO 
AfLtK)  = A(K,K) 

A(K,K>  = T 

TF  (A3S(OEAL(T))  *■  ABF<  AIMflGd)  ) » 15 1 13»  15 

1 7 lP(^n  = *' 

GO  TO  6*’ 

15  00  20  T = KFl,  M 

A<I,K)  = -A(T,K)/T 
OONTIMU^ 

00  4T-  J = KPl,  N 
T = A (L , J) 

A{L,J>  = A<K-,J) 

A(<,J)  = T 

IF  ( ARS(RFAL  (T)  » *■  AnS(  AIMAG  (T»  M 25  t 40»  25 
?5  00  3C  I = KPl,  N 

A(I,J)  = A(I,J)  *■  T*A<I,K) 

3r  nONTTMlJF 

UP  CONTTNUP 

fi'’  oonttnuf 

65  IF  (A(N,N>  .EO.  (0.0,3. 0)1  IPCNl  = G 
RPTMPN 
ENTPy  LEQS3 

SOLUTION  OF  A REAL  LIN-A®  STSTEM  TRI ANGULARIZEO  BY  SUBROUTINE  LUo 

NMl  = N - 1 
IF  (NMl)  2,  9'',  67 

67  00  K = 1,  NMl 

KDi  = «-  ♦ 1 

L = IP(K) 

T = B(L) 

9(LI  = 9(K) 

''(K)  = T 

00  69  T = KPl,  N 

B(T»  = 8(1)  *■  T*A(I,K) 

69  OONTINU'^ 

7P  CONTINUE 

00  A**  L = 1,  NHl 
KMl  = N - L 
K = YMl  ♦ « 

B(K)  = P(K)/A(K,K) 

T = -f’(K) 

no  75  T = 1 , KMl 

n( I ) = 0(1)  + t^a  (I, K) 

75  continue 

S'*  CONTINUE 

9P  P(l)  = n(l>/A (1,1) 

RETIJPN 

ENO 


54 


REFERENCES 


1.  Rice,  E.  J.,  “Attenuation  of  Sound  in  Soft  Walled  Circular  Ducts," 
AFOSR-UTIAS  Symposium  on  Aerodynamic  Noise,  University  of 
Toronto  Press,  1969. 

2.  Baumeister,  K.  J.,  "Application  of  Finite  Difference  Techniques 
to  Noise  Propagation  in  Jet  Engine  Ducts,"  NASA  TM  X-68261 ; also 
presented  at  ASME  Winter  Annual  Meeting,  Detroit,  Mich.,  November 
1973. 

3.  Baumeister,  K.  J., and  Rice,  E.  J.,  "A  Difference  Theory  for  Noise 
Propagation  in  an  Acoustically  Lined  Duct  with  Mean  Flow,"  Progress 
in  Astronautics  and  Aeronautics-Aeroacoustics:  Jet  and  Combustion 

Noise;  Duct  Acoustics,  Vol . 37,  edited  by  H.  T.  Nagamatsu, 

J.  V.  O'Keefe,  and  I.  R,  Schwartz,  AIAA,  New  York,  1975,  pp.  435-453. 

4.  Rice,  E.  J.,  "Propagation  of  Waves  in  an  Acoustically  Lined  Duct  with 
a Mean  Flow,"  Basic  Aerodynamic  Noise  Research,  NASA  SP-207,  1969, 
pp.  345-355. 

5.  Nelsen,  M.  D.,  Linscheid,  L.  L.,  Dinwiddie  III,  B.  A., and  Hall,  Jr., 
0.  J.,  "Study  and  Development  of  Acoustic  Treatment  for  Jet  Engine 
Tailpipes,"  NASA  CR  1853,  1971,  pp.  19-23. 

6.  Quinn,  D.  W.,  "A  Finite  Difference  Method  for  Computing  Sound 
Propagation  in  Nonuniform  Ducts,"  AIAA  Journal , Vol.  13,  September 
1975. 

7.  Isaacson,  E.,and  Keller,  H.  B.,  Analysis  of  Numerical  Methods, 

Wiley,  New  York,  1966,  Chap.  9. 


55 

<tU.S.Gt>vern merit  Printing  OMIce:  1979  — 657-002/86 


